library(tidyverse) # tidy universe
library(kableExtra) # html-table
library(patchwork) # combine plots
library(lmerTest) # mixed model
library(car) # ANOVA
library(emmeans) # post-hoc
library(performance) # model performance
library(Hmisc) # Computes correlation and p-value matrix
library(viridis) # colour scaleCorrelations + Mixed models (Igf & animal welfare)
R Libraries
set.seed(1989)Data
Read data
Data processing done in file “igf-biomarker-summary-stats.qmd”.
load("./data/data-processed.RData")Correlations of parameters
All pairwise correlations
dat.cor.all <- dat.w |>
dplyr::select(where(is.numeric))Calculate correlations and p values:
cor.all <- rcorr(as.matrix(dat.cor.all))
cor.all_matrix <- round(cor.all$r, 2)
pval.all_matrix <- round(cor.all$P, 3)cor.all_matrix |>
kable(caption = "Person correlation coefficients") |>
kable_styling(bootstrap_options = c("striped", "hover"),
font_size = 12) |>
scroll_box()| insem.age | total.piglets | prop.males | prop.stillborn | prop.spreizer | prop.lbw | mean.bodyweight | cort.ser.105dpc | cort.ser.8dpp | cort.sal.105dpc | cort.sal.8dpp | bioact.ser.30dpc | bioact.ser.105dpc | bioact.ser.8dpp | igf2.ser.105dpc | igf2.ser.8dpp | igf1.ser.105dpc | igf1.ser.8dpp | igfbp2.ser.105dpc | igfbp2.ser.8dpp | igfbp3.ser.105dpc | igfbp3.ser.8dpp | proteolysis.105dpc | proteolysis.8dpp | stc1.105dpc | stc1.8dpp | calc.30dpc | calc.105dpc | calc.8dpp | igfbp23.ser.105dpc | igfbp23.ser.8dpp | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| insem.age | 1.00 | 0.35 | -0.06 | 0.25 | -0.21 | -0.20 | 0.16 | -0.19 | -0.44 | 0.02 | -0.19 | -0.29 | -0.20 | -0.21 | 0.07 | -0.28 | 0.53 | 0.42 | 0.18 | 0.32 | -0.15 | -0.10 | 0.19 | 0.13 | 0.15 | 0.10 | -0.10 | 0.00 | 0.01 | 0.12 | 0.05 |
| total.piglets | 0.35 | 1.00 | 0.19 | 0.33 | 0.31 | 0.29 | -0.40 | 0.09 | -0.02 | 0.23 | -0.01 | 0.13 | 0.05 | 0.16 | -0.14 | -0.40 | -0.05 | -0.04 | 0.14 | 0.12 | -0.23 | -0.11 | 0.17 | -0.30 | -0.26 | -0.30 | 0.20 | 0.32 | 0.15 | 0.21 | 0.05 |
| prop.males | -0.06 | 0.19 | 1.00 | 0.09 | 0.07 | 0.21 | -0.21 | 0.24 | 0.49 | -0.09 | -0.11 | -0.03 | -0.12 | -0.16 | 0.20 | 0.18 | -0.16 | -0.21 | 0.04 | 0.18 | 0.25 | 0.29 | 0.08 | -0.04 | -0.23 | -0.22 | 0.06 | 0.07 | -0.05 | -0.16 | -0.12 |
| prop.stillborn | 0.25 | 0.33 | 0.09 | 1.00 | -0.10 | -0.11 | 0.00 | -0.09 | 0.43 | -0.06 | 0.33 | -0.22 | -0.07 | -0.11 | 0.06 | -0.15 | -0.19 | -0.28 | 0.08 | 0.24 | 0.01 | 0.09 | 0.37 | 0.17 | -0.38 | -0.40 | -0.05 | 0.09 | 0.00 | -0.11 | -0.13 |
| prop.spreizer | -0.21 | 0.31 | 0.07 | -0.10 | 1.00 | 0.27 | -0.25 | 0.05 | -0.10 | 0.12 | 0.05 | 0.21 | 0.15 | 0.27 | 0.09 | -0.24 | -0.05 | -0.04 | 0.20 | 0.09 | -0.11 | 0.05 | -0.16 | -0.16 | -0.44 | -0.39 | 0.30 | 0.32 | -0.15 | 0.16 | 0.03 |
| prop.lbw | -0.20 | 0.29 | 0.21 | -0.11 | 0.27 | 1.00 | -0.88 | -0.11 | -0.18 | 0.43 | -0.05 | 0.18 | 0.10 | 0.12 | 0.37 | 0.28 | -0.04 | 0.28 | 0.14 | -0.12 | 0.27 | -0.03 | 0.15 | -0.42 | -0.01 | -0.06 | -0.07 | -0.03 | -0.34 | 0.01 | 0.03 |
| mean.bodyweight | 0.16 | -0.40 | -0.21 | 0.00 | -0.25 | -0.88 | 1.00 | -0.13 | 0.01 | -0.42 | -0.04 | -0.24 | -0.15 | -0.13 | -0.40 | -0.30 | -0.05 | -0.27 | -0.14 | 0.08 | -0.13 | 0.11 | -0.25 | 0.53 | 0.00 | 0.06 | 0.04 | -0.03 | 0.28 | -0.09 | -0.10 |
| cort.ser.105dpc | -0.19 | 0.09 | 0.24 | -0.09 | 0.05 | -0.11 | -0.13 | 1.00 | 0.49 | 0.31 | -0.13 | 0.24 | 0.15 | 0.13 | 0.26 | -0.15 | 0.23 | 0.13 | 0.14 | 0.35 | 0.07 | 0.15 | -0.34 | 0.20 | 0.06 | -0.06 | 0.23 | 0.54 | 0.24 | -0.01 | 0.14 |
| cort.ser.8dpp | -0.44 | -0.02 | 0.49 | 0.43 | -0.10 | -0.18 | 0.01 | 0.49 | 1.00 | -0.22 | 0.02 | -0.40 | -0.44 | -0.53 | 0.11 | 0.09 | -0.09 | -0.16 | 0.11 | 0.40 | 0.35 | 0.23 | 0.02 | 0.33 | -0.08 | -0.10 | -0.25 | 0.06 | -0.08 | -0.24 | -0.18 |
| cort.sal.105dpc | 0.02 | 0.23 | -0.09 | -0.06 | 0.12 | 0.43 | -0.42 | 0.31 | -0.22 | 1.00 | -0.33 | 0.51 | 0.48 | 0.54 | 0.40 | -0.29 | 0.33 | 0.43 | 0.03 | 0.13 | -0.27 | -0.08 | 0.35 | -0.30 | 0.09 | -0.10 | 0.32 | 0.24 | 0.16 | 0.47 | 0.32 |
| cort.sal.8dpp | -0.19 | -0.01 | -0.11 | 0.33 | 0.05 | -0.05 | -0.04 | -0.13 | 0.02 | -0.33 | 1.00 | -0.35 | -0.33 | -0.25 | -0.15 | 0.35 | -0.03 | -0.21 | -0.19 | -0.48 | 0.29 | 0.08 | -0.49 | 0.40 | 0.04 | 0.13 | 0.04 | -0.02 | -0.05 | -0.43 | -0.29 |
| bioact.ser.30dpc | -0.29 | 0.13 | -0.03 | -0.22 | 0.21 | 0.18 | -0.24 | 0.24 | -0.40 | 0.51 | -0.35 | 1.00 | 0.90 | 0.84 | 0.00 | -0.41 | 0.21 | -0.16 | 0.11 | -0.11 | -0.34 | -0.27 | 0.01 | -0.34 | -0.22 | -0.27 | 0.32 | 0.34 | 0.37 | 0.34 | 0.37 |
| bioact.ser.105dpc | -0.20 | 0.05 | -0.12 | -0.07 | 0.15 | 0.10 | -0.15 | 0.15 | -0.44 | 0.48 | -0.33 | 0.90 | 1.00 | 0.89 | 0.13 | -0.46 | 0.31 | -0.07 | 0.20 | -0.09 | -0.29 | -0.25 | 0.07 | -0.18 | -0.34 | -0.38 | 0.31 | 0.32 | 0.26 | 0.39 | 0.36 |
| bioact.ser.8dpp | -0.21 | 0.16 | -0.16 | -0.11 | 0.27 | 0.12 | -0.13 | 0.13 | -0.53 | 0.54 | -0.25 | 0.84 | 0.89 | 1.00 | -0.04 | -0.50 | 0.23 | 0.02 | 0.15 | -0.11 | -0.32 | -0.20 | 0.03 | -0.24 | -0.25 | -0.29 | 0.50 | 0.41 | 0.30 | 0.43 | 0.35 |
| igf2.ser.105dpc | 0.07 | -0.14 | 0.20 | 0.06 | 0.09 | 0.37 | -0.40 | 0.26 | 0.11 | 0.40 | -0.15 | 0.00 | 0.13 | -0.04 | 1.00 | 0.13 | 0.38 | 0.47 | 0.50 | 0.31 | 0.01 | -0.16 | 0.23 | -0.07 | -0.13 | -0.23 | -0.16 | -0.14 | -0.43 | 0.44 | 0.29 |
| igf2.ser.8dpp | -0.28 | -0.40 | 0.18 | -0.15 | -0.24 | 0.28 | -0.30 | -0.15 | 0.09 | -0.29 | 0.35 | -0.41 | -0.46 | -0.50 | 0.13 | 1.00 | -0.10 | 0.03 | -0.11 | -0.15 | 0.35 | 0.27 | 0.00 | -0.24 | 0.53 | 0.64 | -0.21 | -0.21 | -0.17 | -0.41 | -0.34 |
| igf1.ser.105dpc | 0.53 | -0.05 | -0.16 | -0.19 | -0.05 | -0.04 | -0.05 | 0.23 | -0.09 | 0.33 | -0.03 | 0.21 | 0.31 | 0.23 | 0.38 | -0.10 | 1.00 | 0.66 | 0.09 | 0.14 | -0.36 | -0.58 | -0.37 | -0.29 | 0.28 | 0.15 | 0.06 | -0.02 | -0.23 | 0.62 | 0.87 |
| igf1.ser.8dpp | 0.42 | -0.04 | -0.21 | -0.28 | -0.04 | 0.28 | -0.27 | 0.13 | -0.16 | 0.43 | -0.21 | -0.16 | -0.07 | 0.02 | 0.47 | 0.03 | 0.66 | 1.00 | 0.01 | -0.04 | -0.08 | -0.30 | -0.23 | -0.35 | 0.39 | 0.22 | 0.06 | -0.13 | -0.35 | 0.54 | 0.60 |
| igfbp2.ser.105dpc | 0.18 | 0.14 | 0.04 | 0.08 | 0.20 | 0.14 | -0.14 | 0.14 | 0.11 | 0.03 | -0.19 | 0.11 | 0.20 | 0.15 | 0.50 | -0.11 | 0.09 | 0.01 | 1.00 | 0.68 | 0.27 | 0.15 | 0.10 | -0.26 | -0.43 | -0.46 | -0.16 | -0.05 | -0.27 | 0.22 | 0.16 |
| igfbp2.ser.8dpp | 0.32 | 0.12 | 0.18 | 0.24 | 0.09 | -0.12 | 0.08 | 0.35 | 0.40 | 0.13 | -0.48 | -0.11 | -0.09 | -0.11 | 0.31 | -0.15 | 0.14 | -0.04 | 0.68 | 1.00 | 0.05 | 0.22 | 0.38 | -0.23 | -0.33 | -0.36 | -0.21 | -0.07 | -0.30 | 0.13 | 0.13 |
| igfbp3.ser.105dpc | -0.15 | -0.23 | 0.25 | 0.01 | -0.11 | 0.27 | -0.13 | 0.07 | 0.35 | -0.27 | 0.29 | -0.34 | -0.29 | -0.32 | 0.01 | 0.35 | -0.36 | -0.08 | 0.27 | 0.05 | 1.00 | 0.72 | -0.27 | 0.24 | -0.06 | -0.03 | -0.13 | -0.01 | -0.12 | -0.62 | -0.49 |
| igfbp3.ser.8dpp | -0.10 | -0.11 | 0.29 | 0.09 | 0.05 | -0.03 | 0.11 | 0.15 | 0.23 | -0.08 | 0.08 | -0.27 | -0.25 | -0.20 | -0.16 | 0.27 | -0.58 | -0.30 | 0.15 | 0.22 | 0.72 | 1.00 | -0.09 | 0.43 | -0.09 | -0.02 | 0.16 | 0.20 | 0.05 | -0.58 | -0.68 |
| proteolysis.105dpc | 0.19 | 0.17 | 0.08 | 0.37 | -0.16 | 0.15 | -0.25 | -0.34 | 0.02 | 0.35 | -0.49 | 0.01 | 0.07 | 0.03 | 0.23 | 0.00 | -0.37 | -0.23 | 0.10 | 0.38 | -0.27 | -0.09 | 1.00 | -0.30 | -0.25 | -0.28 | -0.31 | -0.02 | -0.08 | 0.16 | -0.20 |
| proteolysis.8dpp | 0.13 | -0.30 | -0.04 | 0.17 | -0.16 | -0.42 | 0.53 | 0.20 | 0.33 | -0.30 | 0.40 | -0.34 | -0.18 | -0.24 | -0.07 | -0.24 | -0.29 | -0.35 | -0.26 | -0.23 | 0.24 | 0.43 | -0.30 | 1.00 | -0.07 | -0.05 | 0.15 | 0.08 | 0.25 | -0.54 | -0.41 |
| stc1.105dpc | 0.15 | -0.26 | -0.23 | -0.38 | -0.44 | -0.01 | 0.00 | 0.06 | -0.08 | 0.09 | 0.04 | -0.22 | -0.34 | -0.25 | -0.13 | 0.53 | 0.28 | 0.39 | -0.43 | -0.33 | -0.06 | -0.09 | -0.25 | -0.07 | 1.00 | 0.95 | -0.01 | -0.16 | 0.14 | -0.12 | 0.11 |
| stc1.8dpp | 0.10 | -0.30 | -0.22 | -0.40 | -0.39 | -0.06 | 0.06 | -0.06 | -0.10 | -0.10 | 0.13 | -0.27 | -0.38 | -0.29 | -0.23 | 0.64 | 0.15 | 0.22 | -0.46 | -0.36 | -0.03 | -0.02 | -0.28 | -0.05 | 0.95 | 1.00 | 0.03 | -0.13 | 0.12 | -0.23 | -0.07 |
| calc.30dpc | -0.10 | 0.20 | 0.06 | -0.05 | 0.30 | -0.07 | 0.04 | 0.23 | -0.25 | 0.32 | 0.04 | 0.32 | 0.31 | 0.50 | -0.16 | -0.21 | 0.06 | 0.06 | -0.16 | -0.21 | -0.13 | 0.16 | -0.31 | 0.15 | -0.01 | 0.03 | 1.00 | 0.67 | 0.42 | 0.20 | 0.04 |
| calc.105dpc | 0.00 | 0.32 | 0.07 | 0.09 | 0.32 | -0.03 | -0.03 | 0.54 | 0.06 | 0.24 | -0.02 | 0.34 | 0.32 | 0.41 | -0.14 | -0.21 | -0.02 | -0.13 | -0.05 | -0.07 | -0.01 | 0.20 | -0.02 | 0.08 | -0.16 | -0.13 | 0.67 | 1.00 | 0.54 | -0.14 | -0.11 |
| calc.8dpp | 0.01 | 0.15 | -0.05 | 0.00 | -0.15 | -0.34 | 0.28 | 0.24 | -0.08 | 0.16 | -0.05 | 0.37 | 0.26 | 0.30 | -0.43 | -0.17 | -0.23 | -0.35 | -0.27 | -0.30 | -0.12 | 0.05 | -0.08 | 0.25 | 0.14 | 0.12 | 0.42 | 0.54 | 1.00 | -0.12 | -0.15 |
| igfbp23.ser.105dpc | 0.12 | 0.21 | -0.16 | -0.11 | 0.16 | 0.01 | -0.09 | -0.01 | -0.24 | 0.47 | -0.43 | 0.34 | 0.39 | 0.43 | 0.44 | -0.41 | 0.62 | 0.54 | 0.22 | 0.13 | -0.62 | -0.58 | 0.16 | -0.54 | -0.12 | -0.23 | 0.20 | -0.14 | -0.12 | 1.00 | 0.76 |
| igfbp23.ser.8dpp | 0.05 | 0.05 | -0.12 | -0.13 | 0.03 | 0.03 | -0.10 | 0.14 | -0.18 | 0.32 | -0.29 | 0.37 | 0.36 | 0.35 | 0.29 | -0.34 | 0.87 | 0.60 | 0.16 | 0.13 | -0.49 | -0.68 | -0.20 | -0.41 | 0.11 | -0.07 | 0.04 | -0.11 | -0.15 | 0.76 | 1.00 |
pval.all_matrix |>
kable(caption = "p values of Person correlation coefficients") |>
kable_styling(bootstrap_options = c("striped", "hover"),
font_size = 12) |>
scroll_box()| insem.age | total.piglets | prop.males | prop.stillborn | prop.spreizer | prop.lbw | mean.bodyweight | cort.ser.105dpc | cort.ser.8dpp | cort.sal.105dpc | cort.sal.8dpp | bioact.ser.30dpc | bioact.ser.105dpc | bioact.ser.8dpp | igf2.ser.105dpc | igf2.ser.8dpp | igf1.ser.105dpc | igf1.ser.8dpp | igfbp2.ser.105dpc | igfbp2.ser.8dpp | igfbp3.ser.105dpc | igfbp3.ser.8dpp | proteolysis.105dpc | proteolysis.8dpp | stc1.105dpc | stc1.8dpp | calc.30dpc | calc.105dpc | calc.8dpp | igfbp23.ser.105dpc | igfbp23.ser.8dpp | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| insem.age | NA | 0.028 | 0.719 | 0.115 | 0.202 | 0.207 | 0.337 | 0.432 | 0.049 | 0.919 | 0.422 | 0.075 | 0.223 | 0.187 | 0.754 | 0.228 | 0.012 | 0.050 | 0.286 | 0.057 | 0.379 | 0.563 | 0.430 | 0.578 | 0.536 | 0.672 | 0.530 | 0.999 | 0.950 | 0.468 | 0.766 |
| total.piglets | 0.028 | NA | 0.230 | 0.040 | 0.054 | 0.071 | 0.010 | 0.720 | 0.920 | 0.337 | 0.956 | 0.442 | 0.776 | 0.326 | 0.545 | 0.083 | 0.830 | 0.846 | 0.399 | 0.470 | 0.178 | 0.512 | 0.490 | 0.193 | 0.282 | 0.209 | 0.221 | 0.046 | 0.344 | 0.214 | 0.764 |
| prop.males | 0.719 | 0.230 | NA | 0.589 | 0.657 | 0.189 | 0.200 | 0.328 | 0.029 | 0.728 | 0.634 | 0.878 | 0.467 | 0.323 | 0.401 | 0.457 | 0.469 | 0.342 | 0.834 | 0.284 | 0.143 | 0.090 | 0.757 | 0.867 | 0.336 | 0.362 | 0.736 | 0.662 | 0.757 | 0.352 | 0.489 |
| prop.stillborn | 0.115 | 0.040 | 0.589 | NA | 0.555 | 0.491 | 0.981 | 0.720 | 0.056 | 0.801 | 0.153 | 0.177 | 0.666 | 0.494 | 0.809 | 0.530 | 0.392 | 0.200 | 0.635 | 0.151 | 0.956 | 0.587 | 0.116 | 0.477 | 0.111 | 0.091 | 0.746 | 0.600 | 0.985 | 0.520 | 0.460 |
| prop.spreizer | 0.202 | 0.054 | 0.657 | 0.555 | NA | 0.091 | 0.115 | 0.825 | 0.687 | 0.624 | 0.848 | 0.209 | 0.355 | 0.092 | 0.705 | 0.311 | 0.840 | 0.873 | 0.254 | 0.588 | 0.538 | 0.758 | 0.511 | 0.496 | 0.060 | 0.103 | 0.068 | 0.046 | 0.347 | 0.360 | 0.852 |
| prop.lbw | 0.207 | 0.071 | 0.189 | 0.491 | 0.091 | NA | 0.000 | 0.664 | 0.436 | 0.070 | 0.842 | 0.281 | 0.527 | 0.473 | 0.110 | 0.235 | 0.857 | 0.214 | 0.421 | 0.477 | 0.117 | 0.857 | 0.550 | 0.066 | 0.961 | 0.820 | 0.653 | 0.843 | 0.030 | 0.958 | 0.857 |
| mean.bodyweight | 0.337 | 0.010 | 0.200 | 0.981 | 0.115 | 0.000 | NA | 0.582 | 0.979 | 0.071 | 0.858 | 0.147 | 0.368 | 0.421 | 0.079 | 0.203 | 0.829 | 0.216 | 0.406 | 0.644 | 0.435 | 0.530 | 0.306 | 0.016 | 0.987 | 0.795 | 0.823 | 0.861 | 0.084 | 0.617 | 0.563 |
| cort.ser.105dpc | 0.432 | 0.720 | 0.328 | 0.720 | 0.825 | 0.664 | 0.582 | NA | 0.032 | 0.208 | 0.604 | 0.323 | 0.540 | 0.582 | 0.300 | 0.565 | 0.354 | 0.590 | 0.556 | 0.137 | 0.785 | 0.527 | 0.339 | 0.565 | 0.826 | 0.813 | 0.353 | 0.016 | 0.326 | 0.974 | 0.567 |
| cort.ser.8dpp | 0.049 | 0.920 | 0.029 | 0.056 | 0.687 | 0.436 | 0.979 | 0.032 | NA | 0.366 | 0.942 | 0.079 | 0.055 | 0.017 | 0.659 | 0.712 | 0.699 | 0.489 | 0.649 | 0.080 | 0.131 | 0.326 | 0.958 | 0.292 | 0.746 | 0.671 | 0.286 | 0.812 | 0.737 | 0.301 | 0.446 |
| cort.sal.105dpc | 0.919 | 0.337 | 0.728 | 0.801 | 0.624 | 0.070 | 0.071 | 0.208 | 0.366 | NA | 0.172 | 0.026 | 0.036 | 0.017 | 0.098 | 0.246 | 0.168 | 0.069 | 0.912 | 0.601 | 0.263 | 0.741 | 0.325 | 0.367 | 0.732 | 0.692 | 0.184 | 0.331 | 0.524 | 0.042 | 0.186 |
| cort.sal.8dpp | 0.422 | 0.956 | 0.634 | 0.153 | 0.848 | 0.842 | 0.858 | 0.604 | 0.942 | 0.172 | NA | 0.132 | 0.152 | 0.287 | 0.552 | 0.139 | 0.912 | 0.376 | 0.413 | 0.032 | 0.221 | 0.737 | 0.130 | 0.200 | 0.874 | 0.583 | 0.882 | 0.947 | 0.841 | 0.058 | 0.216 |
| bioact.ser.30dpc | 0.075 | 0.442 | 0.878 | 0.177 | 0.209 | 0.281 | 0.147 | 0.323 | 0.079 | 0.026 | 0.132 | NA | 0.000 | 0.000 | 1.000 | 0.072 | 0.341 | 0.485 | 0.518 | 0.539 | 0.043 | 0.121 | 0.966 | 0.143 | 0.362 | 0.257 | 0.050 | 0.035 | 0.021 | 0.046 | 0.030 |
| bioact.ser.105dpc | 0.223 | 0.776 | 0.467 | 0.666 | 0.355 | 0.527 | 0.368 | 0.540 | 0.055 | 0.036 | 0.152 | 0.000 | NA | 0.000 | 0.577 | 0.039 | 0.165 | 0.761 | 0.233 | 0.595 | 0.085 | 0.148 | 0.764 | 0.436 | 0.154 | 0.113 | 0.055 | 0.044 | 0.106 | 0.018 | 0.031 |
| bioact.ser.8dpp | 0.187 | 0.326 | 0.323 | 0.494 | 0.092 | 0.473 | 0.421 | 0.582 | 0.017 | 0.017 | 0.287 | 0.000 | 0.000 | NA | 0.882 | 0.024 | 0.293 | 0.920 | 0.384 | 0.511 | 0.057 | 0.232 | 0.901 | 0.303 | 0.300 | 0.231 | 0.001 | 0.009 | 0.063 | 0.009 | 0.036 |
| igf2.ser.105dpc | 0.754 | 0.545 | 0.401 | 0.809 | 0.705 | 0.110 | 0.079 | 0.300 | 0.659 | 0.098 | 0.552 | 1.000 | 0.577 | 0.882 | NA | 0.574 | 0.102 | 0.037 | 0.026 | 0.178 | 0.966 | 0.501 | 0.471 | 0.820 | 0.596 | 0.361 | 0.499 | 0.569 | 0.060 | 0.054 | 0.207 |
| igf2.ser.8dpp | 0.228 | 0.083 | 0.457 | 0.530 | 0.311 | 0.235 | 0.203 | 0.565 | 0.712 | 0.246 | 0.139 | 0.072 | 0.039 | 0.024 | 0.574 | NA | 0.676 | 0.895 | 0.646 | 0.535 | 0.127 | 0.254 | 1.000 | 0.430 | 0.025 | 0.004 | 0.381 | 0.366 | 0.466 | 0.074 | 0.147 |
| igf1.ser.105dpc | 0.012 | 0.830 | 0.469 | 0.392 | 0.840 | 0.857 | 0.829 | 0.354 | 0.699 | 0.168 | 0.912 | 0.341 | 0.165 | 0.293 | 0.102 | 0.676 | NA | 0.001 | 0.677 | 0.530 | 0.105 | 0.005 | 0.211 | 0.310 | 0.253 | 0.540 | 0.806 | 0.924 | 0.297 | 0.002 | 0.000 |
| igf1.ser.8dpp | 0.050 | 0.846 | 0.342 | 0.200 | 0.873 | 0.214 | 0.216 | 0.590 | 0.489 | 0.069 | 0.376 | 0.485 | 0.761 | 0.920 | 0.037 | 0.895 | 0.001 | NA | 0.958 | 0.873 | 0.723 | 0.170 | 0.442 | 0.225 | 0.095 | 0.361 | 0.793 | 0.561 | 0.109 | 0.009 | 0.003 |
| igfbp2.ser.105dpc | 0.286 | 0.399 | 0.834 | 0.635 | 0.254 | 0.421 | 0.406 | 0.556 | 0.649 | 0.912 | 0.413 | 0.518 | 0.233 | 0.384 | 0.026 | 0.646 | 0.677 | 0.958 | NA | 0.000 | 0.114 | 0.367 | 0.699 | 0.282 | 0.066 | 0.049 | 0.373 | 0.785 | 0.113 | 0.188 | 0.341 |
| igfbp2.ser.8dpp | 0.057 | 0.470 | 0.284 | 0.151 | 0.588 | 0.477 | 0.644 | 0.137 | 0.080 | 0.601 | 0.032 | 0.539 | 0.595 | 0.511 | 0.178 | 0.535 | 0.530 | 0.873 | 0.000 | NA | 0.774 | 0.202 | 0.124 | 0.340 | 0.166 | 0.131 | 0.228 | 0.670 | 0.077 | 0.438 | 0.448 |
| igfbp3.ser.105dpc | 0.379 | 0.178 | 0.143 | 0.956 | 0.538 | 0.117 | 0.435 | 0.785 | 0.131 | 0.263 | 0.221 | 0.043 | 0.085 | 0.057 | 0.966 | 0.127 | 0.105 | 0.723 | 0.114 | 0.774 | NA | 0.000 | 0.270 | 0.325 | 0.812 | 0.913 | 0.456 | 0.945 | 0.484 | 0.000 | 0.002 |
| igfbp3.ser.8dpp | 0.563 | 0.512 | 0.090 | 0.587 | 0.758 | 0.857 | 0.530 | 0.527 | 0.326 | 0.741 | 0.737 | 0.121 | 0.148 | 0.232 | 0.501 | 0.254 | 0.005 | 0.170 | 0.367 | 0.202 | 0.000 | NA | 0.711 | 0.064 | 0.715 | 0.920 | 0.359 | 0.251 | 0.782 | 0.000 | 0.000 |
| proteolysis.105dpc | 0.430 | 0.490 | 0.757 | 0.116 | 0.511 | 0.550 | 0.306 | 0.339 | 0.958 | 0.325 | 0.130 | 0.966 | 0.764 | 0.901 | 0.471 | 1.000 | 0.211 | 0.442 | 0.699 | 0.124 | 0.270 | 0.711 | NA | 0.207 | 0.479 | 0.442 | 0.197 | 0.928 | 0.741 | 0.534 | 0.428 |
| proteolysis.8dpp | 0.578 | 0.193 | 0.867 | 0.477 | 0.496 | 0.066 | 0.016 | 0.565 | 0.292 | 0.367 | 0.200 | 0.143 | 0.436 | 0.303 | 0.820 | 0.430 | 0.310 | 0.225 | 0.282 | 0.340 | 0.325 | 0.064 | 0.207 | NA | 0.833 | 0.878 | 0.524 | 0.745 | 0.284 | 0.018 | 0.085 |
| stc1.105dpc | 0.536 | 0.282 | 0.336 | 0.111 | 0.060 | 0.961 | 0.987 | 0.826 | 0.746 | 0.732 | 0.874 | 0.362 | 0.154 | 0.300 | 0.596 | 0.025 | 0.253 | 0.095 | 0.066 | 0.166 | 0.812 | 0.715 | 0.479 | 0.833 | NA | 0.000 | 0.960 | 0.525 | 0.559 | 0.614 | 0.664 |
| stc1.8dpp | 0.672 | 0.209 | 0.362 | 0.091 | 0.103 | 0.820 | 0.795 | 0.813 | 0.671 | 0.692 | 0.583 | 0.257 | 0.113 | 0.231 | 0.361 | 0.004 | 0.540 | 0.361 | 0.049 | 0.131 | 0.913 | 0.920 | 0.442 | 0.878 | 0.000 | NA | 0.909 | 0.594 | 0.616 | 0.341 | 0.791 |
| calc.30dpc | 0.530 | 0.221 | 0.736 | 0.746 | 0.068 | 0.653 | 0.823 | 0.353 | 0.286 | 0.184 | 0.882 | 0.050 | 0.055 | 0.001 | 0.499 | 0.381 | 0.806 | 0.793 | 0.373 | 0.228 | 0.456 | 0.359 | 0.197 | 0.524 | 0.960 | 0.909 | NA | 0.000 | 0.008 | 0.252 | 0.808 |
| calc.105dpc | 0.999 | 0.046 | 0.662 | 0.600 | 0.046 | 0.843 | 0.861 | 0.016 | 0.812 | 0.331 | 0.947 | 0.035 | 0.044 | 0.009 | 0.569 | 0.366 | 0.924 | 0.561 | 0.785 | 0.670 | 0.945 | 0.251 | 0.928 | 0.745 | 0.525 | 0.594 | 0.000 | NA | 0.000 | 0.413 | 0.517 |
| calc.8dpp | 0.950 | 0.344 | 0.757 | 0.985 | 0.347 | 0.030 | 0.084 | 0.326 | 0.737 | 0.524 | 0.841 | 0.021 | 0.106 | 0.063 | 0.060 | 0.466 | 0.297 | 0.109 | 0.113 | 0.077 | 0.484 | 0.782 | 0.741 | 0.284 | 0.559 | 0.616 | 0.008 | 0.000 | NA | 0.472 | 0.396 |
| igfbp23.ser.105dpc | 0.468 | 0.214 | 0.352 | 0.520 | 0.360 | 0.958 | 0.617 | 0.974 | 0.301 | 0.042 | 0.058 | 0.046 | 0.018 | 0.009 | 0.054 | 0.074 | 0.002 | 0.009 | 0.188 | 0.438 | 0.000 | 0.000 | 0.534 | 0.018 | 0.614 | 0.341 | 0.252 | 0.413 | 0.472 | NA | 0.000 |
| igfbp23.ser.8dpp | 0.766 | 0.764 | 0.489 | 0.460 | 0.852 | 0.857 | 0.563 | 0.567 | 0.446 | 0.186 | 0.216 | 0.030 | 0.031 | 0.036 | 0.207 | 0.147 | 0.000 | 0.003 | 0.341 | 0.448 | 0.002 | 0.000 | 0.428 | 0.085 | 0.664 | 0.791 | 0.808 | 0.517 | 0.396 | 0.000 | NA |
Selected correlations
Select columns and merge data frames:
dat.cor.sel <- dat.w |>
dplyr::select(mean.bodyweight, prop.lbw, igf2.ser.105dpc, igf2.ser.8dpp) |>
mutate(prop.lbw = prop.lbw/100) |>
rename(`Piglet birth weight` = mean.bodyweight,
`Proportion LBW` = prop.lbw,
`IGF2 105dpc` = igf2.ser.105dpc,
`IGF2 8dpp` = igf2.ser.8dpp)Calculate correlations and p values:
cor.sel <- rcorr(as.matrix(dat.cor.sel))
cor.sel_matrix <- as_tibble(cor.sel$r, rownames = "Var1")
pval.sel_matrix <- as_tibble(cor.sel$P, rownames = "Var1")cor.sel_matrix |>
kable(caption = "Person correlation coefficients") |>
kable_styling(bootstrap_options = c("striped", "hover"),
font_size = 12) |>
scroll_box()| Var1 | Piglet birth weight | Proportion LBW | IGF2 105dpc | IGF2 8dpp |
|---|---|---|---|---|
| Piglet birth weight | 1.0000000 | -0.8823905 | -0.4013240 | -0.2974563 |
| Proportion LBW | -0.8823905 | 1.0000000 | 0.3686377 | 0.2780085 |
| IGF2 105dpc | -0.4013240 | 0.3686377 | 1.0000000 | 0.1336761 |
| IGF2 8dpp | -0.2974563 | 0.2780085 | 0.1336761 | 1.0000000 |
pval.sel_matrix |>
kable(caption = "p values of Person correlation coefficients") |>
kable_styling(bootstrap_options = c("striped", "hover"),
font_size = 12) |>
scroll_box()| Var1 | Piglet birth weight | Proportion LBW | IGF2 105dpc | IGF2 8dpp |
|---|---|---|---|---|
| Piglet birth weight | NA | 0.0000000 | 0.0794637 | 0.2027779 |
| Proportion LBW | 0.0000000 | NA | 0.1097394 | 0.2353017 |
| IGF2 105dpc | 0.0794637 | 0.1097394 | NA | 0.5742138 |
| IGF2 8dpp | 0.2027779 | 0.2353017 | 0.5742138 | NA |
Heatmap:
# Convert matrices into long format for ggplot
cor_long <- pivot_longer(cor.sel_matrix,
cols = !Var1,
names_to = "Var2",
values_to = "cor.val")
p_long <- pivot_longer(pval.sel_matrix,
cols = !Var1,
names_to = "Var2",
values_to = "p.val")
# Merge correlation values and p-values
heatmap_data <- cor_long |>
left_join(p_long, by = c("Var1", "Var2")) |>
# this makes it to characters
mutate(cor.val = sprintf("%.2f", cor.val),
p.val = sprintf("%.3f", p.val)) |>
mutate(cor.val = case_when(cor.val == "1.00" ~ NA_character_,
TRUE ~ cor.val)) |>
mutate(p.val = case_when(p.val == "0.000" ~ "<0.001",
is.na(p.val) ~ NA_character_,
TRUE ~ paste0("=",p.val))) |>
mutate(label = case_when(is.na(cor.val) ~ NA_character_,
TRUE ~ paste0(cor.val, "\n(p", p.val, ")")))# Create the heatmap
plot.corr <- heatmap_data |>
ggplot(aes(Var1, Var2, fill = as.numeric(cor.val))) +
geom_tile(color = "white") +
scale_fill_viridis(limits = c(-1, 1), guide="none") +
# Overlay correlation & p-values
geom_text(aes(label = label), color = "white", size = 4) +
my_theme_cor +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
labs(title = "",
x = "",
y = "")plot.corr
Statistical model: general design
To correctly account for the nested data structure in the analysis using a linear mixed model with the lmerTest package in R, we need to account for the hierarchical structure in the data:
- Husbandry Level (conventional/ecological): Each sow is housed in either of the two facilities.
- Sow Level: Each sow gave birth 1, 2 or 3 times, so these events are nested within sows.
- Time Point Level: Measurements were taken before and after each birth, introducing repeated measures within each pregnancy.
- control for age of the sows
This type of nested structure is best represented by a model with random effects at the different levels to account for the dependencies within each level. The random effects allow for variation within subjects and
contr = lmerControl(optimizer = "bobyqa",
optCtrl = list(maxfun = 10000000),
calc.derivs = FALSE)Cortisol serum
Model
hist(dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.cort1 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)summary(mod.cort1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.ser")),
value)
Control: contr
REML criterion at convergence: 22.6
Scaled residuals:
Min 1Q Median 3Q Max
-1.45942 -0.49449 -0.08495 0.57811 1.38490
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.02937 0.1714
sow (Intercept) 0.01239 0.1113
Residual 0.05441 0.2333
Number of obs: 39, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 3.40152 0.10154 14.82783 33.498 2.18e-15 ***
husbandryecological -0.03632 0.14635 16.13178 -0.248 0.807
time8dpp -0.10436 0.10432 17.29836 -1.000 0.331
husbandryecological:time8dpp 0.11841 0.15039 17.69062 0.787 0.442
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.694
time8dpp -0.514 0.356
hsbndrycl:8 0.356 -0.532 -0.694round(drop1(mod.cort1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.034 0.034 1 17.691 0.62 0.442Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.cort1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 0.0406 1 0.8404
time 0.3976 1 0.5283
husbandry:time 0.6199 1 0.4311car::Anova(mod.cort1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 0.0388 1 10.877 0.8474
time 0.4022 1 17.439 0.5342
husbandry:time 0.6151 1 17.473 0.4434Model diagnostics
performance::check_model(mod.cort1)
Emmeans & Effect sizes
Most popular effect-size measure is probably Cohen’s d, which is defined as the observed difference, divided by the population SD; and obviously Cohen effect sizes are close cousins of pairwise differences. They are available via the eff_size() function of the emmeans package.
Husbandry
Emmeans:
emm <- emmeans(mod.cort1,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 28.5 2.48 8.61 23.4 34.7
ecological 29.1 2.57 9.27 23.9 35.5
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.977 0.121 8.93 1 -0.185 0.8577
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort1), # Residual standard deviation from model
edf = df.residual(mod.cort1)) # Residual degrees of freedomSince 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.0981 0.532 8.93 -1.3 1.11
Results are averaged over the levels of: time
sigma used for effect sizes: 0.2333
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.cort1,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 29.5 2.16 16.1 25.2 34.4
8dpp 28.2 2.02 14.8 24.2 32.8
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.05 0.0787 17.7 1 0.600 0.5558
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort1),
edf = df.residual(mod.cort1))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.194 0.323 17.7 -0.486 0.874
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.2333
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.cort1,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 30.0 3.05 14.8 24.2 37.3
ecological 28.9 3.05 17.4 23.2 36.1
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 27.0 2.75 14.8 21.8 33.6
ecological 29.3 2.98 14.8 23.6 36.5
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.037 0.152 16.1 1 0.248 0.8071
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.921 0.132 14.8 1 -0.572 0.5762
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort1),
edf = df.residual(mod.cort1))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.156 0.628 16.1 -1.17 1.485
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.352 0.617 14.8 -1.67 0.965
sigma used for effect sizes: 0.2333
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.cort1,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 30.0 3.05 14.8 24.2 37.3
8dpp 27.0 2.75 14.8 21.8 33.6
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 28.9 3.05 17.4 23.2 36.1
8dpp 29.3 2.98 14.8 23.6 36.5
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.110 0.116 17.3 1 1.000 0.3309
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.986 0.107 18.1 1 -0.130 0.8982
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort1),
edf = df.residual(mod.cort1))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.4474 0.451 17.3 -0.502 1.397
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.0602 0.464 18.1 -1.036 0.915
sigma used for effect sizes: 0.2333
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.cort1 <- dat.l |>
dplyr::filter(parameter == "cort.ser") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 60),
breaks = seq(0, 60, 20),
minor_breaks = seq(0, 60, by = 2) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "Cortisol (serum) [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.cort1
Cortisol saliva
Model
hist(dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.cort2 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.cort2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "cort.sal")),
value)
Control: contr
REML criterion at convergence: 101.8
Scaled residuals:
Min 1Q Median 3Q Max
-2.29670 -0.33630 0.06505 0.77409 1.36353
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.000 0.0000
sow (Intercept) 0.000 0.0000
Residual 0.827 0.9094
Number of obs: 39, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.0473 0.3031 35.0000 6.754 7.96e-08 ***
husbandryecological -0.4572 0.4179 35.0000 -1.094 0.2813
time8dpp -0.7537 0.4179 35.0000 -1.804 0.0799 .
husbandryecological:time8dpp 0.9303 0.5831 35.0000 1.595 0.1196
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.725
time8dpp -0.725 0.526
hsbndrycl:8 0.520 -0.717 -0.717
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.cort2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 2.105 2.105 1 35 2.546 0.12Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.cort2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 0.0049 1 0.9439
time 0.8969 1 0.3436
husbandry:time 2.5455 1 0.1106car::Anova(mod.cort2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 0.0058 1 9.7306 0.9409
time 0.8843 1 17.6972 0.3597
husbandry:time 2.5370 1 17.7258 0.1289Model diagnostics
performance::check_model(mod.cort2)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.cort2,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 5.31 1.11 35 3.48 8.12
ecological 5.36 1.09 35 3.55 8.09
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.992 0.289 35 1 -0.027 0.9785
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort2),
edf = df.residual(mod.cort2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.00871 0.321 35 -0.66 0.642
Results are averaged over the levels of: time
sigma used for effect sizes: 0.9094
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.cort2,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 6.16 1.290 35 4.03 9.42
8dpp 4.62 0.939 35 3.06 6.98
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.33 0.389 35 1 0.990 0.3291
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort2),
edf = df.residual(mod.cort2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.317 0.323 35 -0.338 0.973
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.9094
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.cort2,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 7.75 2.35 35 4.19 14.34
ecological 4.90 1.41 35 2.74 8.79
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 3.65 1.05 35 2.03 6.54
ecological 5.85 1.68 35 3.26 10.49
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.580 0.660 35 1 1.094 0.2813
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.623 0.253 35 1 -1.163 0.2526
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort2),
edf = df.residual(mod.cort2))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.503 0.464 35 -0.439 1.444
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.520 0.452 35 -1.438 0.397
sigma used for effect sizes: 0.9094
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.cort2,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 7.75 2.35 35 4.19 14.34
8dpp 3.65 1.05 35 2.03 6.54
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 4.90 1.41 35 2.74 8.79
8dpp 5.85 1.68 35 3.26 10.49
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 2.125 0.888 35 1 1.804 0.0799
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.838 0.341 35 1 -0.434 0.6668
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.cort2),
edf = df.residual(mod.cort2))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.829 0.471 35 -0.127 1.785
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.194 0.448 35 -1.103 0.715
sigma used for effect sizes: 0.9094
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.cort2 <- dat.l |>
dplyr::filter(parameter == "cort.sal") |>
droplevels() |>
drop_na(value) |>
# make plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 30),
breaks = seq(0, 30, 10),
minor_breaks = seq(0, 30, by = 2) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "Cortisol (saliva) [ng/mg protein]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.cort2
IGF bioactivity
Model
hist(dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.bioact <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.bioact)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: droplevels(dplyr::filter(drop_na(dplyr::filter(dat.l, parameter ==
"bioact.ser"), value), time != "30dpc"))
Control: contr
REML criterion at convergence: 143.2
Scaled residuals:
Min 1Q Median 3Q Max
-2.29365 -0.36326 0.07091 0.43812 1.94984
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.3260 0.5709
sow (Intercept) 0.0000 0.0000
Residual 0.1372 0.3703
Number of obs: 80, groups: litter.no:sow, 40; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.11600 0.15217 50.82288 33.620 < 2e-16 ***
husbandryecological -0.76647 0.21520 50.82288 -3.562 0.000811 ***
time8dpp -0.05793 0.11711 38.00000 -0.495 0.623690
husbandryecological:time8dpp 0.20329 0.16562 38.00000 1.227 0.227228
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.707
time8dpp -0.385 0.272
hsbndrycl:8 0.272 -0.385 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.bioact, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.207 0.207 1 38 1.506 0.227Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.bioact,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 11.2029 1 0.0008167 ***
time 0.2786 1 0.5976150
husbandry:time 1.5065 1 0.2196781
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.bioact,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 11.0519 1 11.923 0.00611 **
time 0.2786 1 38.000 0.60068
husbandry:time 1.5065 1 38.000 0.22723
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.bioact)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.bioact,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 161.9 22.7 38 121.8 215
ecological 83.3 11.7 38 62.7 111
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.94 0.386 38 1 3.347 0.0019
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.bioact),
edf = df.residual(mod.bioact))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.8 0.557 38 0.669 2.92
Results are averaged over the levels of: time
sigma used for effect sizes: 0.3703
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.cort2,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 6.16 1.290 35 4.03 9.42
8dpp 4.62 0.939 35 3.06 6.98
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.33 0.389 35 1 0.990 0.3291
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.bioact),
edf = df.residual(mod.bioact))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.779 0.79 35 -0.824 2.38
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.3703
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.bioact,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 166.7 25.4 50.8 122.8 226
ecological 77.4 11.8 50.8 57.1 105
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 157.3 23.9 50.8 115.9 213
ecological 89.6 13.6 50.8 66.0 122
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 2.15 0.463 50.8 1 3.562 0.0008
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.76 0.378 50.8 1 2.617 0.0117
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.bioact),
edf = df.residual(mod.bioact))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 2.07 0.606 50.8 0.853 3.29
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.52 0.595 50.8 0.327 2.71
sigma used for effect sizes: 0.3703
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.bioact,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels(),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 166.7 25.4 50.8 122.8 226
8dpp 157.3 23.9 50.8 115.9 213
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 77.4 11.8 50.8 57.1 105
8dpp 89.6 13.6 50.8 66.0 122
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.060 0.124 38 1 0.495 0.6237
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.865 0.101 38 1 -1.241 0.2222
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.bioact),
edf = df.residual(mod.bioact))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.156 0.316 38 -0.484 0.797
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.392 0.318 38 -1.036 0.251
sigma used for effect sizes: 0.3703
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.bioact <- dat.l |>
dplyr::filter(parameter == "bioact.ser") |>
drop_na(value) |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 600),
breaks = seq(0, 600, 200),
minor_breaks = seq(0, 600, by = 50) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "IGF bioactivity [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.bioact
IGF1 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.igf1 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igf1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf1.ser")),
value)
Control: contr
REML criterion at convergence: 77
Scaled residuals:
Min 1Q Median 3Q Max
-1.90876 -0.52177 -0.03551 0.52631 1.57329
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 2.225e-01 4.717e-01
sow (Intercept) 6.468e-17 8.042e-09
Residual 1.642e-01 4.052e-01
Number of obs: 44, groups: litter.no:sow, 22; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.8938 0.1966 30.0492 24.887 <2e-16 ***
husbandryecological -0.5121 0.2663 30.0492 -1.923 0.0640 .
time8dpp 0.5069 0.1812 20.0000 2.797 0.0111 *
husbandryecological:time8dpp 0.4368 0.2453 20.0000 1.780 0.0902 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.739
time8dpp -0.461 0.340
hsbndrycl:8 0.340 -0.461 -0.739
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igf1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.52 0.52 1 20 3.17 0.09 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.igf1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 1.5445 1 0.21395
time 37.2036 1 1.064e-09 ***
husbandry:time 3.1696 1 0.07502 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igf1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 1.4352 1 11.177 0.25570
time 37.2036 1 20.000 5.826e-06 ***
husbandry:time 3.1696 1 20.000 0.09022 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igf1)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.igf1,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 172 30.0 20 119 247
ecological 128 20.4 20 92 179
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.34 0.317 20 1 1.243 0.2283
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf1),
edf = df.residual(mod.igf1))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.725 0.589 20 -0.504 1.95
Results are averaged over the levels of: time
sigma used for effect sizes: 0.4052
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.igf1,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 103 13.8 30.1 78.7 136
8dpp 213 28.4 30.1 162.6 280
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.484 0.0594 20 1 -5.912 <.0001
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf1),
edf = df.residual(mod.igf1))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.79 0.367 20 -2.56 -1.02
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.4052
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.igf1,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 133 26.2 30.1 89.3 199
ecological 80 14.4 30.1 55.4 115
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 222 43.6 30.1 148.3 331
ecological 205 36.9 30.1 142.4 296
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.67 0.444 30.1 1 1.923 0.0640
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.08 0.287 30.1 1 0.283 0.7793
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf1),
edf = df.residual(mod.igf1))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.264 0.673 30.1 -0.111 2.64
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.186 0.657 30.1 -1.157 1.53
sigma used for effect sizes: 0.4052
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.igf1,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 133 26.2 30.1 89.3 199
8dpp 222 43.6 30.1 148.3 331
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 80 14.4 30.1 55.4 115
8dpp 205 36.9 30.1 142.4 296
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.602 0.1090 20 1 -2.797 0.0111
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.389 0.0644 20 1 -5.705 <.0001
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf1),
edf = df.residual(mod.igf1))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.25 0.47 20 -2.23 -0.27
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -2.33 0.49 20 -3.35 -1.31
sigma used for effect sizes: 0.4052
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.igf1 <- dat.l |>
dplyr::filter(parameter == "igf1.ser") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 500),
breaks = seq(0, 500, 100),
minor_breaks = seq(0, 500, by = 20) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "IGF1 (serum) [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.igf1
IGF2 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.igf2 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igf2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igf2.ser")),
value)
Control: contr
REML criterion at convergence: 15.3
Scaled residuals:
Min 1Q Median 3Q Max
-3.3607 -0.5221 0.0581 0.4752 1.6754
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 1.265e-02 1.125e-01
sow (Intercept) 5.401e-17 7.349e-09
Residual 5.793e-02 2.407e-01
Number of obs: 40, groups: litter.no:sow, 20; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 4.26837 0.08401 34.87974 50.808 < 2e-16 ***
husbandryecological -0.33300 0.11881 34.87974 -2.803 0.008215 **
time8dpp 0.42980 0.10764 18.00000 3.993 0.000853 ***
husbandryecological:time8dpp 0.43848 0.15222 18.00000 2.881 0.009955 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.707
time8dpp -0.641 0.453
hsbndrycl:8 0.453 -0.641 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igf2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.481 0.481 1 18 8.298 0.01 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.igf2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 1.5551 1 0.21239
time 72.7197 1 < 2e-16 ***
husbandry:time 8.2975 1 0.00397 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igf2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 1.4218 1 10.105 0.260359
time 72.7197 1 18.000 9.755e-08 ***
husbandry:time 8.2975 1 18.000 0.009955 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igf2)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.igf2,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 88.5 5.71 18 77.3 101.4
ecological 79.0 5.10 18 69.0 90.5
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.12 0.102 18 1 1.247 0.2284
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf2),
edf = df.residual(mod.igf2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.473 0.383 18 -0.333 1.28
Results are averaged over the levels of: time
sigma used for effect sizes: 0.2407
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.igf2,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 60.5 3.59 34.9 53.6 68.2
8dpp 115.7 6.87 34.9 102.5 130.5
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.523 0.0398 18 1 -8.528 <.0001
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf2),
edf = df.residual(mod.igf2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -2.7 0.458 18 -3.66 -1.73
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.2407
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.igf2,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 71.4 6.00 34.9 60.2 84.7
ecological 51.2 4.30 34.9 43.2 60.7
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 109.7 9.22 34.9 92.5 130.2
ecological 122.0 10.20 34.9 102.8 144.6
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.4 0.166 34.9 1 2.803 0.0082
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.9 0.107 34.9 1 -0.888 0.3807
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf2),
edf = df.residual(mod.igf2))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.384 0.522 34.9 0.323 2.44
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.438 0.497 34.9 -1.446 0.57
sigma used for effect sizes: 0.2407
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.igf2,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 71.4 6.00 34.9 60.2 84.7
8dpp 109.7 9.22 34.9 92.5 130.2
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 51.2 4.30 34.9 43.2 60.7
8dpp 122.0 10.20 34.9 102.8 144.6
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.651 0.0700 18 1 -3.993 0.0009
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.420 0.0452 18 1 -8.067 <.0001
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igf2),
edf = df.residual(mod.igf2))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.79 0.498 18 -2.83 -0.739
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -3.61 0.630 18 -4.93 -2.284
sigma used for effect sizes: 0.2407
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.igf2 <- dat.l |>
dplyr::filter(parameter == "igf2.ser") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 200),
breaks = seq(0, 200, 50),
minor_breaks = seq(0, 200, by = 10) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "IGF2 (serum) [ng/ml]",
fill = "Husbandry") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5))) +
my_theme +
theme(legend.position = c(0.25, 0.85),
legend.box = "horizontal")Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.plot.igf2
Combined plots: Figure 3
# Combine plots with a designated area for the legend
combined <- (plot.cort1 +
plot.cort2 +
plot.bioact +
plot.igf1) +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20),
legend.position = "bottom",
legend.direction = "horizontal")combined
png("./plots/figure3.png",
width = 300, height = 300, units = "mm",
pointsize = 10, res = 600)
combined
dev.off()png
2 Combined plots: Figure 4
# Combine plots with a designated area for the legend
combined <- (plot.igf2 +
plot.corr) +
#plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20))combined
png("./plots/figure4.png",
width = 300, height = 200, units = "mm",
pointsize = 10, res = 600)
combinedWarning: Removed 4 rows containing missing values or values outside the scale range
(`geom_text()`).dev.off()png
2 Litter numbers
plot.bio.li <- dat.w |>
# plot
mutate(jit = jitter(as.numeric(litter.no), 0.3)) |>
ggplot(aes(y = bioact.ser.105dpc)) +
geom_boxplot(aes(x = litter.no,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = litter.no,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 600),
breaks = seq(0, 600, 200),
minor_breaks = seq(0, 600, by = 50) ) +
scale_x_discrete(labels= c("1st", "2nd", "3rd")) +
labs(x = "Litter number",
y = "IGF bioactivity 105dpc [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "bottom") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.lbw.li <- dat.w |>
# plot
mutate(jit = jitter(as.numeric(litter.no), 0.3)) |>
ggplot(aes(y = prop.lbw)) +
geom_boxplot(aes(x = litter.no,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = litter.no,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 100),
breaks = seq(0, 100, 20),
minor_breaks = seq(0, 100, by = 5) ) +
scale_x_discrete(labels= c("1st", "2nd", "3rd")) +
labs(x = "Litter number",
y = "Proportion LBW [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "bottom") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))Combined plots: Figure 5
# Combine plots with a designated area for the legend
combined <- (plot.bio.li +
plot.lbw.li) +
plot_layout(guides = "collect") +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20),
legend.position = "bottom")combined
png("./plots/figure5.png",
width = 300, height = 200, units = "mm",
pointsize = 10, res = 600)
combined
dev.off()png
2 IGFBP2 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.igfbp2 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)summary(mod.igfbp2)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp2.ser")),
value)
Control: contr
REML criterion at convergence: 63.3
Scaled residuals:
Min 1Q Median 3Q Max
-1.80571 -0.56601 0.02105 0.51087 1.84289
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.07463 0.2732
sow (Intercept) 0.04330 0.2081
Residual 0.05864 0.2422
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.51548 0.11794 14.52980 55.243 <2e-16 ***
husbandryecological -0.46275 0.16544 15.40369 -2.797 0.0133 *
time8dpp -0.03649 0.08072 34.00000 -0.452 0.6541
husbandryecological:time8dpp 0.13669 0.11416 34.00000 1.197 0.2394
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.713
time8dpp -0.342 0.244
hsbndrycl:8 0.242 -0.345 -0.707round(drop1(mod.igfbp2, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.084 0.084 1 34 1.434 0.239Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.igfbp2,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 6.4509 1 0.01109 *
time 0.3114 1 0.57680
husbandry:time 1.4337 1 0.23116
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igfbp2,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 6.3339 1 12.02 0.02704 *
time 0.3114 1 34.00 0.58045
husbandry:time 1.4337 1 34.00 0.23945
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igfbp2)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.igfbp2,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 663 73.5 11.4 520 846
ecological 447 48.6 12.7 353 566
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.48 0.23 12 1 2.540 0.0259
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp2),
edf = df.residual(mod.igfbp2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.63 0.657 12 0.198 3.06
Results are averaged over the levels of: time
sigma used for effect sizes: 0.2422
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.igfbp2,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 536 44.3 15.4 450 639
8dpp 553 45.8 15.4 464 660
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.969 0.0553 34 1 -0.558 0.5805
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp2),
edf = df.residual(mod.igfbp2))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.132 0.236 34 -0.611 0.348
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.2422
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.igfbp2,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 676 79.7 14.5 525 869
ecological 425 49.3 16.4 333 544
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 651 76.8 14.5 506 838
ecological 470 54.5 16.4 368 601
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.59 0.263 15.4 1 2.797 0.0133
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.39 0.229 15.4 1 1.971 0.0670
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp2),
edf = df.residual(mod.igfbp2))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.91 0.703 15.4 0.415 3.41
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.35 0.693 15.4 -0.128 2.82
sigma used for effect sizes: 0.2422
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.igfbp2,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 676 79.7 14.5 525 869
8dpp 651 76.8 14.5 506 838
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 425 49.3 16.4 333 544
8dpp 470 54.5 16.4 368 601
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.037 0.0837 34 1 0.452 0.6541
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.905 0.0730 34 1 -1.241 0.2230
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp2),
edf = df.residual(mod.igfbp2))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.151 0.334 34 -0.527 0.829
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.414 0.335 34 -1.095 0.268
sigma used for effect sizes: 0.2422
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.igfbp2 <- dat.l |>
dplyr::filter(parameter == "igfbp2.ser") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 1250),
breaks = seq(0, 1250, 250),
minor_breaks = seq(0, 1250, by = 100) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "IGF BP2 (serum) [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.igfbp2
IGFBP3 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.igfbp3 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igfbp3)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp3.ser")),
value)
Control: contr
REML criterion at convergence: 159.9
Scaled residuals:
Min 1Q Median 3Q Max
-1.78859 -0.50549 0.09757 0.46868 1.52002
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 5.588e-01 7.475e-01
sow (Intercept) 6.984e-16 2.643e-08
Residual 2.038e-01 4.515e-01
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 6.72944 0.20583 44.24476 32.694 < 2e-16 ***
husbandryecological 0.08475 0.29109 44.24476 0.291 0.77229
time8dpp 0.53588 0.15048 34.00000 3.561 0.00112 **
husbandryecological:time8dpp -0.03664 0.21282 34.00000 -0.172 0.86434
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.707
time8dpp -0.366 0.258
hsbndrycl:8 0.258 -0.366 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igfbp3, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.006 0.006 1 34 0.03 0.864Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.igfbp3,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 0.0601 1 0.8063
time 23.6578 1 1.151e-06 ***
husbandry:time 0.0296 1 0.8633
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igfbp3,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 0.0583 1 10.995 0.8136
time 23.6578 1 34.000 2.577e-05 ***
husbandry:time 0.0296 1 34.000 0.8643
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igfbp3)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.igfbp3,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 1094 210 34 741 1614
ecological 1169 224 34 792 1725
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.936 0.254 34 1 -0.245 0.8078
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp3),
edf = df.residual(mod.igfbp3))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.147 0.6 34 -1.37 1.07
Results are averaged over the levels of: time
sigma used for effect sizes: 0.4515
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.igfbp3,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 873 127 44.2 651 1170
8dpp 1465 213 44.2 1092 1964
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.596 0.0634 34 1 -4.864 <.0001
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp3),
edf = df.residual(mod.igfbp3))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.15 0.256 34 -1.67 -0.626
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.4515
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.igfbp3,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 837 172 44.2 553 1267
ecological 911 187 44.2 602 1379
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 1430 294 44.2 944 2165
ecological 1500 309 44.2 991 2271
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.919 0.267 44.2 1 -0.291 0.7723
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.953 0.277 44.2 1 -0.165 0.8695
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp3),
edf = df.residual(mod.igfbp3))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.188 0.645 44.2 -1.49 1.11
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.107 0.645 44.2 -1.41 1.19
sigma used for effect sizes: 0.4515
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.igfbp3,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 837 172 44.2 553 1267
8dpp 1430 294 44.2 944 2165
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 911 187 44.2 602 1379
8dpp 1500 309 44.2 991 2271
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.585 0.0881 34 1 -3.561 0.0011
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.607 0.0913 34 1 -3.318 0.0022
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp3),
edf = df.residual(mod.igfbp3))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.19 0.349 34 -1.90 -0.477
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -1.11 0.347 34 -1.81 -0.400
sigma used for effect sizes: 0.4515
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.igfbp3 <- dat.l |>
dplyr::filter(parameter == "igfbp3.ser") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 6500),
breaks = seq(0, 6500, 2000),
minor_breaks = seq(0, 6500, by = 200) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "IGF BP3 (serum) [ng/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.igfbp3
IGFBP2/IGFBP3 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.igfbp23 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.igfbp23)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "igfbp23.ser")),
value)
Control: contr
REML criterion at convergence: 144.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.97572 -0.51849 -0.00712 0.49478 1.73489
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.4515 0.6719
sow (Intercept) 0.0000 0.0000
Residual 0.1608 0.4010
Number of obs: 72, groups: litter.no:sow, 36; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) -0.2861 0.1844 44.0488 -1.551 0.127980
husbandryecological -0.5229 0.2608 44.0488 -2.005 0.051155 .
time8dpp -0.5724 0.1337 34.0000 -4.282 0.000143 ***
husbandryecological:time8dpp 0.1733 0.1890 34.0000 0.917 0.365645
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.707
time8dpp -0.362 0.256
hsbndrycl:8 0.256 -0.362 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.igfbp23, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.135 0.135 1 34 0.841 0.366Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.igfbp23,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 3.2202 1 0.07273 .
time 26.4099 1 2.761e-07 ***
husbandry:time 0.8407 1 0.35919
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.igfbp23,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 3.1232 1 10.995 0.1049
time 26.4099 1 34.000 1.134e-05 ***
husbandry:time 0.8407 1 34.000 0.3656
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.igfbp23)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.igfbp23,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 0.564 0.0970 34 0.398 0.800
ecological 0.365 0.0627 34 0.257 0.517
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.55 0.376 34 1 1.795 0.0816
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp23),
edf = df.residual(mod.igfbp23))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.09 0.614 34 -0.159 2.34
Results are averaged over the levels of: time
sigma used for effect sizes: 0.401
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.igfbp23,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 0.578 0.0754 44 0.445 0.752
8dpp 0.356 0.0464 44 0.274 0.463
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.63 0.154 34 1 5.139 <.0001
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp23),
edf = df.residual(mod.igfbp23))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 1.21 0.259 34 0.686 1.74
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.401
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.igfbp23,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 0.751 0.1390 44 0.518 1.089
ecological 0.445 0.0821 44 0.307 0.646
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 0.424 0.0782 44 0.292 0.615
ecological 0.299 0.0551 44 0.206 0.433
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.69 0.44 44 1 2.005 0.0512
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.42 0.37 44 1 1.340 0.1870
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp23),
edf = df.residual(mod.igfbp23))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 1.304 0.660 44 -0.0269 2.64
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.872 0.655 44 -0.4481 2.19
sigma used for effect sizes: 0.401
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.igfbp23,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 0.751 0.1390 44 0.518 1.089
8dpp 0.424 0.0782 44 0.292 0.615
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 0.445 0.0821 44 0.307 0.646
8dpp 0.299 0.0551 44 0.206 0.433
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.77 0.237 34 1 4.282 0.0001
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.49 0.199 34 1 2.986 0.0052
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.igfbp23),
edf = df.residual(mod.igfbp23))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 1.427 0.356 34 0.704 2.15
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.995 0.345 34 0.295 1.70
sigma used for effect sizes: 0.401
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.igfbp23 <- dat.l |>
dplyr::filter(parameter == "igfbp23.ser") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = (value))) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 2.5),
breaks = seq(0, 2.5, 0.5),
minor_breaks = seq(0, 2.5, by = 0.1) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "Molar ratio of IGFBP2/IGFBP3 (serum)",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.igfbp23
Proteolytic activity (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.prot <- lmerTest::lmer(log(value+1) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.prot)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula:
log(value + 1) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "proteolysis")),
value)
Control: contr
REML criterion at convergence: 120.5
Scaled residuals:
Min 1Q Median 3Q Max
-1.4098 -1.1650 0.3079 0.7632 1.3081
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 0.000 0.000
sow (Intercept) 0.000 0.000
Residual 1.413 1.189
Number of obs: 39, groups: litter.no:sow, 20; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 1.6496 0.3963 35.0000 4.163 0.000194 ***
husbandryecological -0.2007 0.5462 35.0000 -0.367 0.715542
time8dpp -0.2647 0.5462 35.0000 -0.485 0.630926
husbandryecological:time8dpp 0.4917 0.7622 35.0000 0.645 0.523064
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.725
time8dpp -0.725 0.526
hsbndrycl:8 0.520 -0.717 -0.717
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.prot, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value + 1) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.588 0.588 1 35 0.416 0.523Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.prot,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value + 1)
Chisq Df Pr(>Chisq)
husbandry 0.0185 1 0.8918
time 0.0010 1 0.9744
husbandry:time 0.4162 1 0.5189car::Anova(mod.prot,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value + 1)
F Df Df.res Pr(>F)
husbandry 0.0175 1 10.932 0.8972
time 0.0007 1 17.691 0.9796
husbandry:time 0.4121 1 17.726 0.5291Model diagnostics
performance::check_model(mod.prot)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.prot,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 3.56 1.25 35 1.62 6.94
ecological 3.77 1.27 35 1.78 7.18
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log(mu + 1) scale
$contrasts
contrast estimate SE df t.ratio p.value
conventional - ecological -0.0452 0.381 35 -0.119 0.9063
Results are averaged over the levels of: time
Note: contrasts are still on the log(mu + 1) scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: satterthwaite Effect sizes:
eff_size(emm,
sigma = sigma(mod.prot),
edf = df.residual(mod.prot))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.038 0.321 35 -0.689 0.613
Results are averaged over the levels of: time
sigma used for effect sizes: 1.189
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.prot,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 3.71 1.29 35 1.70 7.20
8dpp 3.62 1.23 35 1.69 6.93
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log(mu + 1) scale
$contrasts
contrast estimate SE df t.ratio p.value
105dpc - 8dpp 0.0189 0.381 35 0.050 0.9608
Results are averaged over the levels of: husbandry
Note: contrasts are still on the log(mu + 1) scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: satterthwaite Effect sizes:
eff_size(emm,
sigma = sigma(mod.prot),
edf = df.residual(mod.prot))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.0159 0.321 35 -0.635 0.667
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 1.189
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.prot,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 4.21 2.06 35 1.328 10.64
ecological 3.26 1.60 35 0.985 8.14
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 2.99 1.50 35 0.862 7.57
ecological 4.34 2.01 35 1.491 10.46
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log(mu + 1) scale
$contrasts
time = 105dpc:
contrast estimate SE df t.ratio p.value
conventional - ecological 0.201 0.546 35 0.367 0.7155
time = 8dpp:
contrast estimate SE df t.ratio p.value
conventional - ecological -0.291 0.532 35 -0.547 0.5875
Note: contrasts are still on the log(mu + 1) scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: satterthwaite Effect sizes:
eff_size(emm,
sigma = sigma(mod.prot),
edf = df.residual(mod.prot))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.169 0.460 35 -0.765 1.103
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.245 0.448 35 -1.155 0.665
sigma used for effect sizes: 1.189
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.prot,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 4.21 2.06 35 1.328 10.64
8dpp 2.99 1.50 35 0.862 7.57
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 3.26 1.60 35 0.985 8.14
8dpp 4.34 2.01 35 1.491 10.46
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log(mu + 1) scale
$contrasts
husbandry = conventional:
contrast estimate SE df t.ratio p.value
105dpc - 8dpp 0.265 0.546 35 0.485 0.6309
husbandry = ecological:
contrast estimate SE df t.ratio p.value
105dpc - 8dpp -0.227 0.532 35 -0.427 0.6720
Note: contrasts are still on the log(mu + 1) scale. Consider using
regrid() if you want contrasts of back-transformed estimates.
Degrees-of-freedom method: satterthwaite Effect sizes:
eff_size(emm,
sigma = sigma(mod.prot),
edf = df.residual(mod.prot))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.223 0.460 35 -0.712 1.157
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.191 0.448 35 -1.100 0.718
sigma used for effect sizes: 1.189
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.prot <- dat.l |>
dplyr::filter(parameter == "proteolysis") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 21),
breaks = seq(0, 21, 5),
minor_breaks = seq(0, 21, by = 1) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "Proteolytic activity (serum) [%]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.prot
Stanniocalcin, STC1 (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.stc1 <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.stc1)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dat.l, parameter == "stc1")),
value)
Control: contr
REML criterion at convergence: 88.2
Scaled residuals:
Min 1Q Median 3Q Max
-1.58308 -0.48078 0.00956 0.47032 1.39341
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 3.17571 1.7821
sow (Intercept) 0.00000 0.0000
Residual 0.05628 0.2372
Number of obs: 38, groups: litter.no:sow, 19; sow, 14
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 5.62613 0.59926 17.29863 9.388 3.32e-08 ***
husbandryecological 1.99293 0.82602 17.29863 2.413 0.0272 *
time8dpp -0.31643 0.11184 17.00000 -2.829 0.0116 *
husbandryecological:time8dpp 0.08471 0.15416 17.00000 0.550 0.5898
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.725
time8dpp -0.093 0.068
hsbndrycl:8 0.068 -0.093 -0.725
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.stc1, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.017 0.017 1 17 0.302 0.59Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.stc1,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 6.1244 1 0.0133326 *
time 12.4733 1 0.0004128 ***
husbandry:time 0.3020 1 0.5826404
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1car::Anova(mod.stc1,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 5.4881 1 10.123 0.040863 *
time 12.4733 1 17.000 0.002561 **
husbandry:time 0.3020 1 17.000 0.589787
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Model diagnostics
performance::check_model(mod.stc1)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.stc1,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 237 141 17 67.3 834
ecological 1814 1030 17 549.5 5987
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.131 0.107 17 1 -2.475 0.0242
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.stc1),
edf = df.residual(mod.stc1))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -8.58 3.63 17 -16.2 -0.912
Results are averaged over the levels of: time
sigma used for effect sizes: 0.2372
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.stc1,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 752 311 17.3 315 1795
8dpp 572 236 17.3 239 1365
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.32 0.101 17 1 3.556 0.0024
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.stc1),
edf = df.residual(mod.stc1))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 1.16 0.356 17 0.403 1.91
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.2372
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.stc1,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 278 166 17.3 78.5 981
ecological 2037 1160 17.3 614.7 6747
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 202 121 17.3 57.2 715
ecological 1615 918 17.3 487.6 5352
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.136 0.113 17.3 1 -2.413 0.0272
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.125 0.103 17.3 1 -2.515 0.0220
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.stc1),
edf = df.residual(mod.stc1))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -8.40 3.64 17.3 -16.1 -0.727
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -8.76 3.66 17.3 -16.5 -1.056
sigma used for effect sizes: 0.2372
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.stc1,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 278 166 17.3 78.5 981
8dpp 202 121 17.3 57.2 715
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 2037 1160 17.3 614.7 6747
8dpp 1615 918 17.3 487.6 5352
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.37 0.153 17 1 2.829 0.0116
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.26 0.134 17 1 2.184 0.0433
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.stc1),
edf = df.residual(mod.stc1))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 1.334 0.501 17 0.27695 2.39
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.977 0.464 17 -0.00245 1.96
sigma used for effect sizes: 0.2372
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.stc1 <- dat.l |>
dplyr::filter(parameter == "stc1") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 10000),
breaks = seq(0, 10000, 2000),
minor_breaks = seq(0, 10000, by = 400) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "STC1 (serum) [pg/ml]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.stc1
Combined plots: Figure 6
# Combine plots with a designated area for the legend
combined <- (plot.igfbp23 +
plot.stc1) +
plot_annotation(tag_levels = 'A') &
theme(plot.tag = element_text(face = "bold", size = 20),
legend.position = "bottom")combined
png("./plots/figure6.png",
width = 300, height = 200, units = "mm",
pointsize = 10, res = 600)
combined
dev.off()png
2 Calcium (serum)
Model
hist(dat.l |>
dplyr::filter(parameter == "calc") |>
droplevels() |>
drop_na(value) |>
pull("value"),
breaks = 30)
hist(log(dat.l |>
dplyr::filter(parameter == "calc") |>
droplevels() |>
drop_na(value) |>
pull("value")),
breaks = 30)
Singularity issues if nested structure defined like (1 | sow/litter.no).
mod.calc <- lmerTest::lmer(log(value) ~
husbandry +
time +
# interaction term
husbandry : time +
# random intercept for sows and for each litter within sow
(1 | sow/litter.no),
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
REML = TRUE,
control = contr)boundary (singular) fit: see help('isSingular')summary(mod.calc)Linear mixed model fit by REML. t-tests use Satterthwaite's method [
lmerModLmerTest]
Formula: log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Data: drop_na(droplevels(dplyr::filter(dplyr::filter(dat.l, parameter ==
"calc"), time != "30dpc")), value)
Control: contr
REML criterion at convergence: 4.1
Scaled residuals:
Min 1Q Median 3Q Max
-5.9530 -0.3184 0.0786 0.3944 1.6808
Random effects:
Groups Name Variance Std.Dev.
litter.no:sow (Intercept) 6.379e-18 2.526e-09
sow (Intercept) 3.982e-02 1.995e-01
Residual 3.860e-02 1.965e-01
Number of obs: 80, groups: litter.no:sow, 40; sow, 15
Fixed effects:
Estimate Std. Error df t value Pr(>|t|)
(Intercept) 2.42516 0.08737 15.46761 27.756 1.29e-14 ***
husbandryecological -0.08834 0.12130 16.18553 -0.728 0.477
time8dpp -0.02237 0.06213 62.42643 -0.360 0.720
husbandryecological:time8dpp 0.10889 0.08786 62.42643 1.239 0.220
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Correlation of Fixed Effects:
(Intr) hsbndr tm8dpp
hsbndryclgc -0.720
time8dpp -0.356 0.256
hsbndrycl:8 0.251 -0.362 -0.707
optimizer (bobyqa) convergence code: 0 (OK)
boundary (singular) fit: see help('isSingular')round(drop1(mod.calc, test = 'Chisq'), 3)Single term deletions using Satterthwaite's method:
Model:
log(value) ~ husbandry + time + husbandry:time + (1 | sow/litter.no)
Sum Sq Mean Sq NumDF DenDF F value Pr(>F)
husbandry:time 0.059 0.059 1 62.426 1.536 0.22Without significant interactions, choose type = 2. With significant interactions, choose type = 3.
car::Anova(mod.calc,
test.statistic = "Chisq",
type = 2)Analysis of Deviance Table (Type II Wald chisquare tests)
Response: log(value)
Chisq Df Pr(>Chisq)
husbandry 0.0899 1 0.7644
time 0.5333 1 0.4652
husbandry:time 1.5360 1 0.2152car::Anova(mod.calc,
test.statistic = "F",
type = 2)Analysis of Deviance Table (Type II Wald F tests with Kenward-Roger df)
Response: log(value)
F Df Df.res Pr(>F)
husbandry 0.0897 1 12.849 0.7694
time 0.5333 1 38.000 0.4697
husbandry:time 1.5360 1 38.000 0.2228Model diagnostics
performance::check_model(mod.calc)
Emmeans & Effect sizes
Husbandry
Emmeans:
emm <- emmeans(mod.calc,
pairwise ~ husbandry,
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
husbandry response SE df lower.CL upper.CL
conventional 11.2 0.913 11.8 9.35 13.4
ecological 10.8 0.845 12.7 9.12 12.8
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.03 0.117 12.2 1 0.300 0.7694
Results are averaged over the levels of: time
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.calc),
edf = df.residual(mod.calc))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.173 0.576 12.2 -1.08 1.42
Results are averaged over the levels of: time
sigma used for effect sizes: 0.1965
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point
Emmeans:
emm <- emmeans(mod.calc,
pairwise ~ time,
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")NOTE: Results may be misleading due to involvement in interactionsemm$emmeans
time response SE df lower.CL upper.CL
105dpc 10.8 0.656 16.2 9.51 12.3
8dpp 11.2 0.677 16.2 9.82 12.7
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.968 0.0425 62.4 1 -0.730 0.4680
Results are averaged over the levels of: husbandry
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.calc),
edf = df.residual(mod.calc))Since 'object' is a list, we are using the contrasts already present. contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.163 0.224 62.4 -0.611 0.284
Results are averaged over the levels of: husbandry
sigma used for effect sizes: 0.1965
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Husbandry:Time point
Emmeans:
emm <- emmeans(mod.calc,
pairwise ~ husbandry|time,
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
time = 105dpc:
husbandry response SE df lower.CL upper.CL
conventional 11.3 0.988 15.5 9.39 13.6
ecological 10.3 0.871 17.0 8.67 12.4
time = 8dpp:
husbandry response SE df lower.CL upper.CL
conventional 11.1 0.966 15.5 9.18 13.3
ecological 11.3 0.949 17.0 9.45 13.5
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
time = 105dpc:
contrast ratio SE df null t.ratio p.value
conventional / ecological 1.09 0.133 16.2 1 0.728 0.4768
time = 8dpp:
contrast ratio SE df null t.ratio p.value
conventional / ecological 0.98 0.119 16.2 1 -0.169 0.8675
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.calc),
edf = df.residual(mod.calc))Since 'object' is a list, we are using the contrasts already present.time = 105dpc:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) 0.450 0.619 16.2 -0.86 1.76
time = 8dpp:
contrast estimate SE df lower.CL upper.CL
(conventional - ecological) -0.105 0.617 16.2 -1.41 1.20
sigma used for effect sizes: 0.1965
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Time point:Husbandry
Emmeans:
emm <- emmeans(mod.calc,
pairwise ~ time|husbandry,
data = dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value),
adjust = "tukey",
lmer.df = "satterthwaite",
type = "response")
emm$emmeans
husbandry = conventional:
time response SE df lower.CL upper.CL
105dpc 11.3 0.988 15.5 9.39 13.6
8dpp 11.1 0.966 15.5 9.18 13.3
husbandry = ecological:
time response SE df lower.CL upper.CL
105dpc 10.3 0.871 17.0 8.67 12.4
8dpp 11.3 0.949 17.0 9.45 13.5
Degrees-of-freedom method: satterthwaite
Confidence level used: 0.95
Intervals are back-transformed from the log scale
$contrasts
husbandry = conventional:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 1.023 0.0635 62.4 1 0.360 0.7201
husbandry = ecological:
contrast ratio SE df null t.ratio p.value
105dpc / 8dpp 0.917 0.0570 62.4 1 -1.393 0.1686
Degrees-of-freedom method: satterthwaite
Tests are performed on the log scale Effect sizes:
eff_size(emm,
sigma = sigma(mod.calc),
edf = df.residual(mod.calc))Since 'object' is a list, we are using the contrasts already present.husbandry = conventional:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) 0.114 0.316 62.4 -0.518 0.746
husbandry = ecological:
contrast estimate SE df lower.CL upper.CL
(105dpc - 8dpp) -0.440 0.318 62.4 -1.077 0.196
sigma used for effect sizes: 0.1965
Degrees-of-freedom method: inherited from satterthwaite when re-gridding
Confidence level used: 0.95 Plot
plot.calc <- dat.l |>
dplyr::filter(parameter == "calc") |>
dplyr::filter(time != "30dpc") |>
droplevels() |>
drop_na(value) |>
# plot
mutate(jit = jitter(as.numeric(time), 0.3)) |>
ggplot(aes(y = value)) +
geom_boxplot(aes(x = time,
fill = husbandry),
col = "black",
outlier.shape = NA, width = 0.5) +
geom_jitter(aes(x = time,
fill = husbandry),
col = "black",
position = position_jitterdodge(jitter.width = 0.15,
dodge.width = 0.5),
size = 2, alpha = 0.7,
show.legend = FALSE) +
scale_fill_manual(labels = c("Conventional",
"Ecological"),
values = c("#FFA040", "#008000")) +
scale_y_continuous(lim = c(0, 21),
breaks = seq(0, 20, 5),
minor_breaks = seq(0, 20, by = 2) ) +
scale_x_discrete(labels= c("105dpc", "8dpp")) +
labs(x = "",
y = "Calcium (serum) [mg/dl]",
fill = "Husbandry") +
my_theme +
theme(legend.position = "top") +
guides(y = guide_axis(minor.ticks = TRUE),
fill = guide_legend(override.aes = list(linetype = 0, size = 5)))plot.calc
How to cite R
“All analyses were performed using R Statistical Software (version 4.4.2; R Core Team 2024)”.
Reference: R Core Team (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. URL https://www.R-project.org/.
citation()To cite R in publications use:
R Core Team (2024). _R: A Language and Environment for Statistical
Computing_. R Foundation for Statistical Computing, Vienna, Austria.
<https://www.R-project.org/>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {R: A Language and Environment for Statistical Computing},
author = {{R Core Team}},
organization = {R Foundation for Statistical Computing},
address = {Vienna, Austria},
year = {2024},
url = {https://www.R-project.org/},
}
We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.version$version.string[1] "R version 4.4.2 (2024-10-31 ucrt)"citation("tidyverse")Um Paket 'tidyverse' in Publikationen zu zitieren, nutzen Sie bitte:
Wickham H, Averick M, Bryan J, Chang W, McGowan LD, François R,
Grolemund G, Hayes A, Henry L, Hester J, Kuhn M, Pedersen TL, Miller
E, Bache SM, Müller K, Ooms J, Robinson D, Seidel DP, Spinu V,
Takahashi K, Vaughan D, Wilke C, Woo K, Yutani H (2019). "Welcome to
the tidyverse." _Journal of Open Source Software_, *4*(43), 1686.
doi:10.21105/joss.01686 <https://doi.org/10.21105/joss.01686>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {Welcome to the {tidyverse}},
author = {Hadley Wickham and Mara Averick and Jennifer Bryan and Winston Chang and Lucy D'Agostino McGowan and Romain François and Garrett Grolemund and Alex Hayes and Lionel Henry and Jim Hester and Max Kuhn and Thomas Lin Pedersen and Evan Miller and Stephan Milton Bache and Kirill Müller and Jeroen Ooms and David Robinson and Dana Paige Seidel and Vitalie Spinu and Kohske Takahashi and Davis Vaughan and Claus Wilke and Kara Woo and Hiroaki Yutani},
year = {2019},
journal = {Journal of Open Source Software},
volume = {4},
number = {43},
pages = {1686},
doi = {10.21105/joss.01686},
}citation("kableExtra")Um Paket 'kableExtra' in Publikationen zu zitieren, nutzen Sie bitte:
Zhu H (2024). _kableExtra: Construct Complex Table with 'kable' and
Pipe Syntax_. R package version 1.4.0,
<https://CRAN.R-project.org/package=kableExtra>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {kableExtra: Construct Complex Table with 'kable' and Pipe Syntax},
author = {Hao Zhu},
year = {2024},
note = {R package version 1.4.0},
url = {https://CRAN.R-project.org/package=kableExtra},
}citation("patchwork")Um Paket 'patchwork' in Publikationen zu zitieren, nutzen Sie bitte:
Pedersen T (2024). _patchwork: The Composer of Plots_. R package
version 1.3.0, <https://CRAN.R-project.org/package=patchwork>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {patchwork: The Composer of Plots},
author = {Thomas Lin Pedersen},
year = {2024},
note = {R package version 1.3.0},
url = {https://CRAN.R-project.org/package=patchwork},
}citation("lmerTest")To cite lmerTest in publications use:
Kuznetsova A, Brockhoff PB, Christensen RHB (2017). "lmerTest
Package: Tests in Linear Mixed Effects Models." _Journal of
Statistical Software_, *82*(13), 1-26. doi:10.18637/jss.v082.i13
<https://doi.org/10.18637/jss.v082.i13>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {{lmerTest} Package: Tests in Linear Mixed Effects Models},
author = {Alexandra Kuznetsova and Per B. Brockhoff and Rune H. B. Christensen},
journal = {Journal of Statistical Software},
year = {2017},
volume = {82},
number = {13},
pages = {1--26},
doi = {10.18637/jss.v082.i13},
}citation("car")To cite the car package in publications use:
Fox J, Weisberg S (2019). _An R Companion to Applied Regression_,
Third edition. Sage, Thousand Oaks CA.
<https://www.john-fox.ca/Companion/>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Book{,
title = {An {R} Companion to Applied Regression},
edition = {Third},
author = {John Fox and Sanford Weisberg},
year = {2019},
publisher = {Sage},
address = {Thousand Oaks {CA}},
url = {https://www.john-fox.ca/Companion/},
}citation("emmeans")Um Paket 'emmeans' in Publikationen zu zitieren, nutzen Sie bitte:
Lenth R (2024). _emmeans: Estimated Marginal Means, aka Least-Squares
Means_. R package version 1.10.6,
<https://CRAN.R-project.org/package=emmeans>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {emmeans: Estimated Marginal Means, aka Least-Squares Means},
author = {Russell V. Lenth},
year = {2024},
note = {R package version 1.10.6},
url = {https://CRAN.R-project.org/package=emmeans},
}citation("performance")Um Paket 'performance' in Publikationen zu zitieren, nutzen Sie bitte:
Lüdecke et al., (2021). performance: An R Package for Assessment,
Comparison and Testing of Statistical Models. Journal of Open Source
Software, 6(60), 3139. https://doi.org/10.21105/joss.03139
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Article{,
title = {{performance}: An {R} Package for Assessment, Comparison and Testing of Statistical Models},
author = {Daniel Lüdecke and Mattan S. Ben-Shachar and Indrajeet Patil and Philip Waggoner and Dominique Makowski},
year = {2021},
journal = {Journal of Open Source Software},
volume = {6},
number = {60},
pages = {3139},
doi = {10.21105/joss.03139},
}citation("Hmisc")Um Paket 'Hmisc' in Publikationen zu zitieren, nutzen Sie bitte:
Harrell Jr F (2025). _Hmisc: Harrell Miscellaneous_. R package
version 5.2-2, <https://CRAN.R-project.org/package=Hmisc>.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {Hmisc: Harrell Miscellaneous},
author = {Frank E {Harrell Jr}},
year = {2025},
note = {R package version 5.2-2},
url = {https://CRAN.R-project.org/package=Hmisc},
}citation("viridis")To cite viridis/viridisLite in publications use:
Simon Garnier, Noam Ross, Robert Rudis, Antônio P. Camargo, Marco
Sciaini, and Cédric Scherer (2024). viridis(Lite) -
Colorblind-Friendly Color Maps for R. viridis package version 0.6.5.
Ein BibTeX-Eintrag für LaTeX-Benutzer ist
@Manual{,
title = {{viridis(Lite)} - Colorblind-Friendly Color Maps for R},
author = {{Garnier} and {Simon} and {Ross} and {Noam} and {Rudis} and {Robert} and {Camargo} and Antônio Pedro and {Sciaini} and {Marco} and {Scherer} and {Cédric}},
year = {2024},
note = {viridis package version 0.6.5},
url = {https://sjmgarnier.github.io/viridis/},
doi = {10.5281/zenodo.4679423},
}Session Info
sessionInfo()R version 4.4.2 (2024-10-31 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)
Matrix products: default
locale:
[1] LC_COLLATE=German_Germany.utf8 LC_CTYPE=German_Germany.utf8
[3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
[5] LC_TIME=German_Germany.utf8
time zone: Europe/Berlin
tzcode source: internal
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] viridis_0.6.5 viridisLite_0.4.2 Hmisc_5.2-2 performance_0.13.0
[5] emmeans_1.10.6 car_3.1-3 carData_3.0-5 lmerTest_3.1-3
[9] lme4_1.1-36 Matrix_1.7-1 patchwork_1.3.0 kableExtra_1.4.0
[13] lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4
[17] purrr_1.0.2 readr_2.1.5 tidyr_1.3.1 tibble_3.2.1
[21] ggplot2_3.5.1 tidyverse_2.0.0
loaded via a namespace (and not attached):
[1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.2.0
[4] TH.data_1.1-3 bayestestR_0.15.1 rpart_4.1.23
[7] digest_0.6.37 timechange_0.3.0 estimability_1.5.1
[10] lifecycle_1.0.4 cluster_2.1.6 survival_3.7-0
[13] magrittr_2.0.3 compiler_4.4.2 rlang_1.1.4
[16] tools_4.4.2 yaml_2.3.10 data.table_1.16.4
[19] knitr_1.49 labeling_0.4.3 htmlwidgets_1.6.4
[22] xml2_1.3.6 abind_1.4-8 multcomp_1.4-26
[25] withr_3.0.2 foreign_0.8-87 numDeriv_2016.8-1.1
[28] datawizard_1.0.0 nnet_7.3-19 grid_4.4.2
[31] xtable_1.8-4 colorspace_2.1-1 scales_1.3.0
[34] MASS_7.3-61 insight_1.0.1 cli_3.6.3
[37] mvtnorm_1.3-3 rmarkdown_2.29 reformulas_0.4.0
[40] generics_0.1.3 rstudioapi_0.17.1 tzdb_0.4.0
[43] minqa_1.2.8 splines_4.4.2 parallel_4.4.2
[46] base64enc_0.1-3 vctrs_0.6.5 boot_1.3-31
[49] sandwich_3.1-1 jsonlite_1.8.9 hms_1.1.3
[52] ggrepel_0.9.6 pbkrtest_0.5.3 htmlTable_2.4.3
[55] Formula_1.2-5 systemfonts_1.1.0 see_0.9.0
[58] glue_1.8.0 nloptr_2.1.1 codetools_0.2-20
[61] stringi_1.8.4 gtable_0.3.6 munsell_0.5.1
[64] pillar_1.10.1 htmltools_0.5.8.1 R6_2.6.1
[67] Rdpack_2.6.2 evaluate_1.0.3 lattice_0.22-6
[70] backports_1.5.0 rbibutils_2.3 broom_1.0.7
[73] Rcpp_1.0.13-1 checkmate_2.3.2 gridExtra_2.3
[76] svglite_2.1.3 coda_0.19-4.1 nlme_3.1-166
[79] mgcv_1.9-1 xfun_0.50 zoo_1.8-12
[82] pkgconfig_2.0.3